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In  a front- end  detection  system,  when  a strong  source 
(earthquakes  as  specifically  discussed  here)  generates  a great  amount  of 
energy  which  lasts  a substantial  length  of  time,  a detector  can  be  jammed 
and  loses  its  detection  capability  or  yields  excessive  false  alarms.  This 
problem  is  more  pronounced  in  a siesmic  detection  system  as  earthquakes 
occur  within  a distance  of  2000  km  from  the  sensors.  Our  attempts  to 
solve  the  problem  are  being  made  by  employing  a number  of  techniques  ap- 
plicable as  a front-end  detector  or  as  a post-detection  processor,  respective- 
ly. 

By  employing  array  beamforming  techniques  which  we  develop- 
ed, we  will  evaluate  the  feasibility  of  using  arrays  to  separate  mixed  signals 
which  are  weak  signals  occurring  after  strong  signals.  They  are  difficult  to 
detect  due  to  the  coda  or  scattered  energy  of  the  preceding  stronger  signal 
interfering  or  overshadowing  the  later  arriving  weak  signal.  The  problem 
here  is  how  to  reject  the  high -amplitude  interfering  wavetrains  of  the  stronger 
earthquake  and  how  to  extract  the  weaker  hidden  signal  from  the  interfering 
environment. 

Two  array  beamforming  techniques  were  studied  for  this  pur- 
pose: one  is  in  the  time  domain  and  can  be  used  as  a front-end  detector  and  the 
other  in  the  frequency  domain  and  intended  for  a post-detection  processor. 

First,  we  applied  to  the  problem  the  Lj  norm  adaptive  beamformer  (ABF)  which 
is  of  time-varying  adaptation  rate  in  terms  of  norm  formulation.  The  ABF 
is  a time-domain  processor  and  can  operate  in  a real-time  and  front-end  detec- 
tion system.  Secondly,  using  the  formulation  of  the  variational  principle  we 
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tried  to  null-steer  the  array  by  the  method  of  weighting  the  array  sensors  so 
that  it  could  reject  the  stronger  interfering  wavetrains  and  extract  the  hidden 
weak  signal  through  an  array  beamforming  process.  This  formulation  was  in 
the  frequency  domain  and  the  technique  can  be  used  for  a post-detection  pro- 
cessor. 

The  ABF  algorithm  described  by  Shen  (1977)  was  modified  in 
this  report.  Mathematical  algorithms  for  the  array- shading  technique  were 
formulated,  software  programs  were  developed  for  the  computations,  and 
experimental  tests  on  recorded  data  from  the  Korean  Seismic  Research  Sta- 
tion (KSRS)  short-period  array  were  conducted.  Section  II  presents  a brief 
description  of  the  adaptive  beamforming  algorithm.  Section  III  describes  the 
physical  formulation  of  array  weighting  algorithms  and  their  computational 
results.  Experimental  procedures  and  results  from  the  recorded-data  tests 
are  presented  in  Section  IV,  followed  by  conclusions  and  suggestions  for 
further  study  in  Section  V.  Finally,  references  are  listed  in  Section  VI. 
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SECTION  II 

THE  ADAPTIVE  BEAMFORMEK  — FRONT-END  DETECTOR 
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A time-domain  linear  adaptive  beamformer  (ABF)  was  develop- 
ed for  use  in  an  operational  front-end  detection  system.  The  objective  is  to 
increase  the  detection  capability  of  a siesmic  array  and  to  extract  the  low- 
magnitude  signals  of  unknown  waveforms  from  noise.  The  algorithm  mini- 
mizes the  squared  filter  output  subject  to  the  constraints  which  pass  energy 
from  steered  location.  Mathematically,  let  y(t)  be  the  adaptive  filter  output 
at  time  t for  an  array  of  N sensors  (channels)  with  a 2J+1  points,  or 
tappings  of  filter  length  per  channel.  Then,  we  write 


y(t)  = 


N J 

£ £ 

j = l n = -J 


a.(n)  X.(t-nAt) 
J J 


(II- 1) 


where  is  time-aligned  data  input  for  the  j 


th 


channel,  At  is  the  sampling 

" th 

interval,  and  a.(n)  is  the  adaptive  filter  coefficient  for  the  j channel 


th  3 

and  n tapping  (filter  point)  of  the  sensor.  In  the  algorithm,  a.(n)  is 

J 


iterated,  or  updated  in  the  time-domain  by  the  expression 


a. 


i‘+At(n)  = a*(n)  + fi\(t)y(t)  f X(t-nAt)  - X.(t-nAt)"!  (II-2) 

J J J J 


where  a^(n)  i*  initialized  to  satisfy  the  Levin's  constraints  (Levin,  1964) 


N 


£ ai 


(n)  = <J 


no 


(II- 3) 


j=l 


II- 1 


i 


for  a center  point  output  described  in  equation  (II-l)  and  where 


Y(t-nAt)  = 


77  2 xj<t-n4t» 

j=l 


(n-4) 


X is  introduced  to  preserve  the  constraints  which  were  initialized  at  the  very 
beginning  of  processing.  The  algorithm  was  formulated  (with  A(t)  = 1. 0)  by 
Burget  al.  (1967)  and  has  been  investigated  by  various  authors  (Kobayashi, 
1970;  Frost,  1972;  Owsley,  1973;  and  Cangi  and  Byun,  1976).  The  algorithm 
achieves  a substantial  noise  reduction  gain  relative  to  the  beamsteer  (time- 
align  and  average)  processor.  However,  experimental  work  on  the  recorded 
seismic  data  indicated  a substantial  signal  loss  compared  with  the  beamsteer 
output.  That  resulted  in  the  optimum  signal-to-noise  ratio  (SNR)  gain  around 
1 dB  relative  to  the  beamsteer  processor. 

A better  time-domain  adaptation  is  needed  to  achieve  a better 
SNR  gain.  It  was  found  that  a quasi  norm  adaptive  filter  improved  the 
performance  both  in  reducing  the  noise  and  in  preventing  the  desired  signal 
from  being  degraded  (Shen,  1977).  Currently,  the  Lj  norm  ABF  has  been 
modified  in  an  attempt  to  make  the  ABF  processor  implementable  for  optimum 
operation  in  various  levels  of  noise  environment.  Presently,  the  time-varying 
adaptation  rate,  A(t)  is  chosen  as  follows: 


7 LTS^ 

A(t)  = ySTS(t)  x p(t) 


(II- 5) 


where 


yLTs(t)  = (1-a)yLTs(t-1)  + “ly(t)l 


(H- 6) 


ySTS(t)  = d-^)ySTS(t-i) + ^|y(t)| 


(II- 7) 


IlM 


-v  -r..>  , 


p(t)  = ^ X^t-nAt)  (II-  8) 

j=l  n=-j 


« is  a long-term  smoothing  constant  of  filter  output  for  noise  adaptation  so 
that  yLTS<t)  can  characterize  the  noise  level  in  the  environment.  fi  is  the 
short-term  smoothing  constant  of  filter  output  intended  to  prevent  the  signal 
loss  when  Ygpgft)  increases.  The  introduction  of  an<*  ^STS^ 

is  intended  to  fix  p for  the  detector  so  that  it  can  operate  in  different  noise 
environments.  It  is  noted  that  the  performance  of  the  adaptive  beamformer 
is  highly  dependent  on  <*»  fit  and  ^ . 
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SECTION  m 


AMPLITUDE  SHADING  OF  ARRAY  SENSORS  — 

POST -DETECTION  PROCESSOR 

A.  INTRODUCTORY  REMARKS 

The  Dolph-Chebyshev  technique  is  often  applied  to  compute 
the  coefficients  for  array  sensors.  It  gives  the  minimum  width  of  main- 
lobe  for  a given  response  of  minor  lobes  or  gives  the  minimum  response  of 
minor  lobes  for  a given  width  of  mainlobe.  The  purpose  of  the  computation 
is  to  obtain  an  array  directional  response  so  that  a signal  being  searched 
or  estimated  can  be  enhanced  from  the  noise  environment.  The  technique 
was  first  developed  by  Dolph  (1946)  who  expanded  the  pattern  function  to 
the  appropriate  Chebyshev  polynomials  and  obtained  an  analytical  solution 
of  the  coefficients  (or  weightings)  for  array  sensors.  However,  the  analytical 
approach  has  the  disadvantage  that  it  is  limited  to  simple  configurations  of 
sensors  such  as  equally  spaced  line  arrays,  rectangular  arrays,  or  circular 
arrays. 

For  the  mixed-signal  separation  problem  which  we  addressed 
in  this  report,  we  tried  to  obtain  the  directional  steering  of  an  array  to  the 
desired  signal  and  simultaneously  attempted  to  create  a 'null'  response  for 
the  array  in  the  direction  of  an  unwanted  interfering  signal.  By  doing  so, 
the  desired  signal  could  pass  the  array  beamforming  processor  without 
being  distorted,  and  on  the  other  hand,  the  interfering  energy  could  be 
reduced  or  eliminated  by  filtered  array  beamforming.  Because  it  is  not 
designed  to  void  the  array  response  in  a specific  direction,  the  Dolph- 
Chebyshev  method  can  not  be  applied  in  this  case,  even  if  simple  geometry 
of  the  array  is  employed.  Furthermore,  the  Dolph-Chebyshev  method, 
which  minimizes  the  maximum  side  lobe  response  of  arrays  of  simple 
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geometry,  cannot  be  applied  to  arbitrary  sensor  configurations  which  typify 
seismic  arrays.  In  light  of  the  difficulties  the  analytical  method  encountered, 
one  has  to  rely  on  a more  general  numerical  approach.  A digital  computer 
perfectly  serves  the  purpose  of  determining  the  null-steer  coefficients  re- 
quired to  optimize  mixed -signal  separation. 

In  this  section,  a numerical  method  will  be  demonstrated  in 
which  we  applied  the  variational  principle  to  minimize  the  total  array  power 
response  subject  to  filter  weight  constraints.  The  physical  model  we  used  is 
an  isotropic  noise  and  two  propagating  signals  - one  for  the  signal  being  sought 
and  the  other  for  the  stronger  interfering  signal  to  be  eliminated. 

B.  PHYSICAL  FORMULATION 

1.  Real  Amplitude  Weighting  Algorithm 

a.  Amplitude  and  Power  Pattern  Received  by  an  Array 

The  configuration  of  array  sensors  considered  here  is  two- 
dimensional  and  is  on  the  earth's  surface.  Let  the  j^1  sensor  of  an  array 
be  located  at  (r^,  <t>^ ) in  polar  coordinates  and  a wave  of  single  frequency 
with  wavenumber  k (=2irf/v,  where  f is  frequency  and  v is  apparent  phase 
velocity)  be  incident  on  the  array  from  azimuthal  angle  <t>  . The  amplitude 
pattern  received  by  the  array  of  N sensors  is  expressed  as 

Eifkr,cos(d>-<f>.)  - kr.cos (<f>  -<6.)1 

e L > > > ° “ J (IQ-1) 

j=l 

where 

i = ^-1  , and 

(f>o  = the  steering  azimuth  of  the  array,  or  the  direction  of 
a desired  signal  to  be  searched  or  estimated. 


m-2 


In  the  seismic  community,  the  azimuthal  angle  is  usually  ref- 
erenced to  the  geographical  North  in  clockwise  rotation,  but  equation  (IU-1) 
should  not  be  changed  because  cosine  is  an  even  function.  Also,  equation 
(III- 1 ) is  the  beamsteer  (evenly  weighted)  amplitude  pattern  of  the  array.  The 
real  amplitude  shading  (weighting)  is  meant  to  multiply  each  array  sensor  by 
a real  number  (or  coefficient)  to  achieve  a certain  design  goal.  Hence  we 
re-write  equation  (m-1)  in  the  form  where  the  j**1  sensor  is  to  be  weighted 

by  a real  number  a.  such  that 

J 

i[kr  coe(<t>-4>  )-  kr  cos(<A  - <6. )] 

F(k = 2^a je  L J J J • (in-2) 

j=l 

Consequently,  the  array  beam  power  pattern  is 

P(k  ,<*>,<£  ) = 
o 


or 

p(k,4>.  4>o)  = 

where 

C.  = cos  I kr.cos(4>- <f>.)  - kr  coB(<f> 

13  l 3 3 j ° J 

- kr^  <(>j  ) + kr^  cos(<£q-  <f>^  )^  . (HI-4) 


F(k,<^> , <ft  ) F*(k,4>,<M 
o o 


^ « i[kr  coa(<f>-<f>  )- kr  co8(<f>  -<f>  ) 

2*2*  Y/e 

/= i j=i 

kr^  cos(<#>-  <f>2  ) + kr^  cos(<^q-  )J  (III- 3a) 


N 


L v!ZZv(c/i 

j=i  l > j 


(in- 3b) 


in- 3 


1 

Normally,  the  configuration  of  a seismic  array  is,  in  some  sense,  arbitrary  ^ 

and  their  sensor  locations  are  expressed  in  the  Cartesian  coordinates  in  which 
we  have 


x.  = r.  sin  d>. 

3 3 3 


and 


y.  = r.  cos  <f>.  . 
J J J 


(HI- 5) 


Again,  we  write  C^.  in  equation  (UI--4)  in  Cartesian  coordinates  here 


C..  = cos 

Zj 


^k(Xj-x^  )(sin4>-  sin <f>Q)  + k(y.-y^  )(cos<£-cos  (Ill- 6) 


b.  Noise  and  Signal  Model 


Power  pattern  received  by  the  array  in  equation  (HI- 3)  is  for 
a single  incident  'ray'  of  signal  or  noise.  For  the  problem  we  addressed  in 
this  report,  we  were  essentially  concerned  with  two  signals  - one  is  the  de- 
sired signal  to  be  enhanced  and  the  other  is  the  interfering  signal  to  be  re- 
duced or  eliminated.  In  addition  to  the  mixed  signals,  we  added  an  isotropic 
noise  to  the  environment  for  practical  realization.  Therefore,  we  have  one 
noise  and  two  signals  as  the  working  physical  model. 


The  energy  for  the  isotropic  noise  should  be  the  integration  of 


equation  (III- 3b)  over  the  whole  two-dimensional  space.  Let  P^(k)  be  the 


power  for  the  noise  component.  Then  we  have 

Zn 


pN(k*<A0)  = f P(k.4>,0o)d4>  . 


(m-7a) 
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Equation  (UI-7a)  can  be  analytically  accomplished  for  certain  configurations 
of  array  geometry  and  was  elegantly  done  for  a line  array  and  a circular  ar- 
ray by  Wang  (1975,  1977).  However,  for  a seismic  array  as  expressed  in 
equation  (III- 3b)  and  (III-6),  we  have  to  rely  on  numerical  integration. 

Practical  applications  and  experiences  suggested  that  energy 

I 

for  the  desired  signal  and  the  interfering  signal  has  to  be  integrated  over  the 
directional  sectors,  respectively,  instead  of  delta  functions  for  single  rays. 
Hence,  we  write 


Pr(k*  <t>) 

L O 


f 


<b  +6 

" x I 


</>  -S 
I I 


P(k A , <f>Q)  d 


(in- 7b) 


where  is  the  azimuth  of  energy  interfering  with  the  signal,  and,  8^  is 
the  integration  range  for  interfering  energy,  and 


ps(k,  *oi  = 


/ 


Vss 


Vss 


p(k,4>,4>o)  d<t> 


(in- 7c) 


for  the  power  of  the  signal  being  enhanced  where  <p  was  defined  before  as  the 

azimuth  for  the  desired  signal  and  8_  is  its  integration  range.  Finally,  the 

o 

total  energy  E(k,  </>q)  being  received  by  the  array  is 


E(k,  </>o)  = PN(k,  <f>Q)  + ggPs(k,  <t>Q)  + gIPI(k>o)  (in- 8) 

where  g and  g are  the  positive  scaling  parameters  for  the  two  mixed  sig- 
b I 

nals.  We  note  that  there  is  no  distribution  function  (in  terms  of  azimuths) 
applied  in  the  integrations  of  equations  (in- 7).  We  assume  that  signal  inten- 
sity is  equally  distributed  over  the  azimuthal  space  for  the  integration  range 
of  interest. 
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c.  Minimization  of  Array  Power  Response 

In  order  to  design  coefficients  which  enhance  the  desired  signal 

and  reduces  the  interfering  energy,  the  array  power  response,  E'(k,  ) must 

o 

be  written  as 


E'(k,  <f>Q)  = PN(k,  4>q)  - ggPs(k,  <t>Q ) + gjP^k,  <l>o)  . (HI- 9) 


We  used  the  variational  principle  with  respect  to  the  parameters  a^,  where 
j = 1,  2,  ...»  N,  so  that  we  minimize  the  power  response  such  that 


dE'(k,<*>o) 


j ~ i # 2,  • . • , N 


(ni-io) 


Recall  that  E'(k,  4>q)  is  a function  of  array  coefficients  a as  indicated  in 
equation  (HI- 3b).  Equation  (III -10)  has  a trivial  solution  which  is  of  no  interest, 
To  avoid  it,  we  imposed  a constraint  condition  on  the  filter  weights 


r *,  ■ 1 • 

j=i 


(in-  ii) 


The  constraints  in  equation  (IH-11)  also  preserve  the  filter  weights  to  have  a 
unity  response  in  array  beam  direction  . 

The  solution  for  equation  (III -10),  along  with  the  constraints  of 
equation  (III-ll)  can  be  obtained  by  Lagrange's  multiplier  method  in  which  we 
write 


L = E'(k ,*Q)+  *(]£  a - 1)  . 


(in- 12) 


in- 6 


Hence,  the  variational  principle  becomes 


and 


dL 

da. 

3 


dL 

d £ 


= 0, 


j = 1 , 2,  ...»  N 


(HI- 31a) 


0 . 


(HI- 13b) 


The  set  of  equations  (HI-1 3a)  (HI-1  3b)  is  an  N+l  simultaneous  equation  and 
their  matrix  form  is 


R1  1 R1  2 
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1 3 


R2  1 R2  2 R2  3 


R3  1 R3  2 


RN  1 RN  2* 


3 3 


V1  N 


2 N 
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3 N 


‘N  N 
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(III- 14) 


where 


Vss 


V*1 


>/  J ■ / c/  jd*-8s  / C/  jd*+  «,  / c/ 

‘V8! 


5) 


t-ss 


C.  . is  referred  to  equation  (HI-6)  and  we  note:  C,  . = 1 when  L - j . 

I J I J 

2.  Complex  Amplitude  Weighting  Algorithm 

In  a straightforward  manner,  equation  (UI-1)  can  be  multiplied 
by  a complex  number  to  achieve  the  shading  purpose  for  both  amplitudes  and 


IH- 7 


■j  ■ 


phases.  Let's  write 


R . IM 
a.  = a.  + la. 

J 3 3 


(III- 16) 


where  j = 1,  2,  ....  N,  for  the  filter  weights.  We  immediately  realized  that 
IM 

a^  was  zero  and  all  filter  weights  were  zero-phased  in  the  case  of  a real  am- 
plitude shading  algorithm.  The  power  for  a single  ray  of  waves  received  by 
the  array  is 


N N 


E Ey<? 

l=l  j=i 


I kfx.-x.Hsin^-sin^  )+  kfy.-y,  )(cos^-cos<^  ) 

L J / ° 3 £ oJ(IH-l7) 


P(k,4>>o)  = (a*)2  + (a* 

J J 

3=1 


IMj2 


- \ 'V  'f.  R R IM  II 

22^L|Vj+a/  ai 

j >i 


) cos(ARG^  ^ 


. IM  R R IM  . .1 

- (a<  a.  - a<a.  ).m(ARG/.)J 


(in-18) 


where  the  argument  ARG  is 


ARG^j  = k(Xj-x^)(sin^-sin^o)  + k(y^-y^  )(cos^-cos^o)  . (Ill- 19) 

We  would  use  the  same  physical  model  described  before  for  the 
formulation  and  accordingly  would  obtain  the  same  generalized  form  of  equa- 
tion (HI-9)  for  the  array  power  response.  Following  the  same  variational 
principle  of  minimizing  the  power  response  with  respect  to  filter  weights. 
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we  required  that 
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and 


«E'(k,  <f>  ) 
o 

« R 
da. 

J 


= 0 


dE'(k,  <f)Q) 

da. 

J 


= 0 


(ni-20) 


(in- 2i) 


where  j = 1,  2,  . . . , N.  To  avoid  the  trivial  solution  where  a.  =0  and 

IM  ^ 

a.  = 0,  it  is  required  to  impose  the  constraints  on  the  filter  weights.  With- 
out getting  into  the  trouble  of  solving  the  second  order  simultaneous  equations, 
the  following  simple  constraints  are  introduced: 


I 


N 


(III- 22) 


(III- 23) 


We  will  show  in  Subsection  III-C  that  the  constraints  in  equations  (III- 22)  and 
(III- 23)  yield  different  solutions  of  filter  weights,  but  result  in  the  same 
power  response  pattern. 
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We  would  use  the  Lagrange's  multiplier  method  again  for  ob- 
taining the  optimum  solution.  Let's  write  the  Lagrangian,  0(k ) in  this  case 


such  that 


0(k)  = E’(k 


for  the  constraints  of  equation  (III-22),  and 


(III- 24) 


0'(k)  = E'(k 


(III- 25) 


for  the  constraints  of  equation  (III- 23).  Then  the  variational  principle  states 


d»(k) 
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(IU-26a) 


(III-  26b) 
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(III-26d) 


where  j = 1,  2,  ...»  N and  where  N is  the  number  of  sensors  in  the  array. 
There  are  2N+2  simultaneous  equations  to  solve  for  2N+2  unknown  variables. 
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The  matrix  form  for  equation  (IH-26)  is  written  here: 
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where  ARC.,  was  defined  in  equation  (III- 19).  We  note  that  ARG  . = 0 if 
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/ = j,  T = T , and  S = -S  . For  the  constraints  of  equation  (HI-23), 
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the  right-hand  side  of  matrix  equation  (HI- 27)  would  be 
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C.  THEORETICALLY  CALCULATED  ARRAY  RESPONSE  PATTERNS 

Computation  for  both  the  real  and  complex  amplitude  weighting 
algorithms  was  done  on  an  IBM  360/44  computer  for  various  physical  design 
parameters  defined  before.  Presented  here  are  only  the  related  examples 
which  were  applied  in  the  test  using  the  recorded  data  from  the  KSRS  short- 
period  seismic  array.  Table  IH-1  shows  the  relative  locations  of  seismic 
sensors  in  the  array. 

Figure  IH-1  shows  the  computed  power  response  pattern  for 
the  beamsteer  (time -alignment  and  average)  processor.  The  array  was 
steered  at  the  290°  azimuth,  the  wave  frequency  was  2.  0 Hz;  and  the  phase 
velocity  was  14.  3 km/ sec.  Velocity  and  azimuth  of  a seismic  bodywave  as 
viewed  by  sensors  in  a given  array  indicated  the  location  of  a seismic  source 
in  the  earth.  In  this  case,  the  corresponding  signal  source  was  in  Tadzhik, 
USSR. 


Figure  HI- 2 shows  the  computed  power  response  pattern  for 
the  real-amplitude-weighted  array.  The  input  parameters  are: 
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TABLE  III-l 

RELATIVE  LOCATIONS  OF  KOREAN  SEISMIC  RESEARCH 
STATIONS  SHORT-PERIOD  ARRAY  SENSORS 
(POSITIVE  Y-AXIS  IS  POINTING  TO  THE  NORTH, 
WHILE  POSITIVE  X TO  THE  EAST) 
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Sensors 

X 

(km) 

Y 

(km) 
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0.  0 

0.  0 
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0.  79 

2.  11 
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2.  08 

0.  13 
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2.  29 

-2.  25 
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-2.  13 
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-4.  81 
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-0.  91 

-3.81 
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-2.  71 
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-0.  95 

18 
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1.  09 
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FIGURE  HI-1 

COMPUTED  BEAMSTEER  RES] — 
2.0  Hz  AND  14.  3 km 
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FIGURE  III- 2 

COMPUTED  POWER  RESPONSE  PATTERN  FOR  REAL- AMPLITUDE 
SHADED  ARRAY  (f  = 2.  0 Hz,  v = 14.  3 km/sec,  g = 0.  0,  g = 0.  0) 


<t>  = 290° 

H - o-o 

gj  = 0.  0 . 
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The  positive  scaling  parameters  for  two  mixed  signals  g and  g were  set 

u I 

to  g = 0.  0 and  g = 0.  0,  suggesting  that  only  the  isotropic  noise  was  used 
b X 

in  the  physical  model.  In  this  figure,  it  was  the  filter  weight  constraints  that 
forced  the  mainlobe  narrower  and  sharper  as  compared  with  the  beamsteer 
response  pattern  in  Figure  III- 1 . The  computed  filter  weights  for  each  sen- 
sor are  also  included  in  Figure  III-2.  This  response  pattern  provides  a test 
for  the  meaning  of  the  constraints  imposed  on  the  filter  weights  (i.  e.  , equa- 
tion (in-11)). 

The  filter  coefficients  in  the  last  figure  were  not  used  in  the 
actual  work  for  solving  the  mixed  signal  separation  problem.  For  that  pur- 
pose, Figure  III- 3 shows  the  computed  power  response  pattern  from  the  filter 
coefficients,  included  in  the  figure,  for  the  suppression  of  the  interfering  sig- 
nal from  the  Ceram  Sea.  The  design  parameters  are: 

4>o  = 290° 

= 182° 
f = 2.  0 Hz 

v = 14.  3 km/sec 


g = 1.0 

gj  = 100.0 

6S  = 15° 


In  this  case,  the  182°  null  steering  azimuth  was  intended  to  reduce  the  inter- 
fering signal  from  a strong  earthquake  in  Ceram,  Indonesia. 

A complex  amplitude  weighting  algorithm  was  also  programmed 
for  computation  and  tested  with  the  recorded  data.  Figure  HI- 4 shows  the 
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FIGURE  HI- 3 

COMPUTED  POWER  RESPONSE  PATTERN  FOR  REAL- AMPLITUDE 
SHADED  ARRAY  FOR  INTERFERING  SIGNAL  AT  182°  AZIMUTH 
(f  = 2.  0 Hz,  v = 14.  3 km/iec,  gg=  1.  0,  gj=  100.  0,  8~=  15°, 

AND  8.=  15°) 


FIGURE  UI-4 


COMPUTED  POWER  RESPONSE  PATTERN  FOR  COMPLEX 
FILTER  WEIGHTS  IN  TABLE  III- 2 (f  = 2.  0 Hz, 
v = 14.  3 km/sec,  g = 0.  0,  AND  g = 0.  0) 


computed  power  response  pattern  for  the  physical  model  of  isotropic  noise 
only.  The  parameters  used  are: 

f = 2.  0 Hz 

v = 14.  3 km/sec 

<b  = 290° 

gs  = 0.0 

gj  = 0.  0 . 
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The  corresponding  computed  filter  weights  were  listed  in  Table  III-2  where 
(a)  was  for  the  constraints  in  equation  (III- 22)  and  (b)  for  the  constraints  in 
equation  (III- 23).  The  two  constraint  conditions  yielded  the  same  response 
pattern  in  Figure  HI-4  but  different  filter  weights.  The  test  in  this  case  pro- 
vided no  actual  application,  but  a practical  realization  for  the  filter  weight 
constraints  used  in  the  algorithm. 

To  demonstrate  the  null  steering  capability  for  the  algorithm. 
Figure  HI- 5 shows  the  computed  power  response  pattern  for  the  complex- 
shaded  array.  The  corresponding  computed  complex  filter  weights  for  the 
two  constraint  conditions  were  listed  in  Table  HI- 3.  The  input  parameters 
for  the  computation  are: 

f = 2.  0 Hz 

v = 14.  3 km/sec 

<t>  = 290° 

o 

= 182° 

*S  = 1<0 

gT  = 125.0 


The  two  constraint  conditions  yielded  the  different  sets  of  filter  weights,  but 
resulted  in  the  same  power  response  pattern  as  mentioned  earlier. 
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COMPLEX  FILTER  WEIGHTS  FOR  FIGURE  III-4 
(f=2.  0 Hz  , v=14.  3 km/sec,  g =0.  0,  g =0.  0) 
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DESIGNED  COMPLEX  FILTER  WEIGHTS  FOR  INTERFERING 
SIGNAL  REDUCTION  AT  182°  AZIMUTH  (f=2.  0 Hz,  v=14.  3 km/sec 
<f>  =290  , <#>=182°,  g =1.0,  g =125.0,  8 =9°,  AND  8 =12°) 
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D.  PRELIMINARY  COMMENTS 

The  power  response  patterns  presented  so  far  were  limited  to 
the  case  of  one  frequency  only.  Limited  computational  experience  on  various 
sets  of  parameters  indicated  that  the  real  amplitude  shading  algorithm  was 
more  stable  than  the  complex  amplitude  shading  algorithm.  The  reason  is 
simple:  The  real -amplitude  algorithm  requires  only  N+l  simultaneous 
equations  for  an  N-channel  array,  while  the  complex-amplitude  algorithm 
needs  2N+ 2 simultaneous  equations.  The  larger  matrix  inversion  yields 
the  poorer  accuracy  for  numerical  solutions,  and  resulted  in  the  less  stable 
system  for  the  algorithm.  Also,  the  present  complex  shading  algorithm 
needed  four  times  as  much  computation  time  as  the  real  amplitude  shading 
algorithm.  Examination  of  the  null-steering  produced  in  the  azimuth  range 
between  160°  and  210°  was  examined  by  comparing  the  maximum  lobes  with- 
in that  range.  For  real  coefficients  the  rejection  in  that  azimuth  range  was 
-40  dB;  for  complex  weights,  -50  dB.  This  compares  to  -16  dB  and  -24  dB, 
respectively,  for  isotropic  noise  rejection  using  real  and  complex  weights. 

Normally,  the  solution  of  filter  weights,  real  or  complex,  de- 
pended on  the  parameters  used;  i.  e.  , f,  v,  g , g , 5 , and  8 . The  use  of 

bib  1 

g functioned  as  a directivity  design  parameter,  while  g as  the  null- 
b X 

steering  parameter  which  controls  the  null- response  levels.  The  stability 
of  the  solution  was  highly  dependent  on  the  azimuthal  separation  between  the 
interfering  signal  and  the  desired  signal,  </>q  - <f> ^ , and  the  wavelength  of  a 
seismic  P wave,  X . 
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EXPERIMENTAL  RESULTS  FROM  RECORDED  SEISMIC  ARRAY  DATA 
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A.  DATA 

The  algorithms  were  tested  using  the  recorded  data  from  the 
short-period  array  at  the  Korean  Seismic  Research  Station  (KSRS).  The 
KSRS  short-period  array  has  19  sensors  with  an  aperture  of  about  10  km. 

The  configuration  consists  of  two  concentric  rings:  6 sensors  in  the  inner 
ring  and  12  sensors  in  the  outer  ring  with  one  in  the  center.  The  relative 
coordinates  of  sensors  have  been  presented  in  Table  III- 1 . The  data  used  in 
this  report  were  digitized  at  sampling  time  interval  At  = 0.  1 seconds. 

Two  earthquake  events  were  used  for  mixed- signal  data:  one 

was  from  Ceram  Sea,  Indonesia  (NEIS:  01/01/77,  19.01.40,  2.  532S,  126.  582 

E,  m^)=  6.  0)  and  the  other  from  Tadzhik,  USSR  (NEIS:  03/08/77,  15.  28.  47, 

38.  03N,  69.  44E,  m,  = 5.  2).  The  first  event  had  a very  strong  and  long  wave- 

fa 

train  recorded  at  KSRS  and  was  used  as  the  strong  interfering  signal.  The 
second  event  had  signal  amplitude  at  the  noise  level  and  was  undetectable  in 
the  single  sensor  data,  but  detectable  in  the  array  beam. 

To  demonstrate  the  performance  of  techniques  we  developed 
for  the  mixed- signal  situations,  the  raw  data  of  the  weaker  signal  were  added 
to  those  of  the  much  stronger  signal.  Figure  IV- 1 shows  the  composite  data 
(in  0.  5 - 3.  5 Hz  passband)  which  simulated  the  recorded  mixed  signal  situa- 
tions. In  this  case,  the  raw  data  of  the  stronger  signal  had  been  scaled  by  a 

factor  of  0.  1 before  they  were  added  to  the  data  of  the  weaker  signal.  The  > 

signal  wavetrains  from  Ceram  Sea  dominate  the  display,  while  the  later  ar- 
riving  signal  from  Tadzhik  is  hidden  by  the  visible  wavetrains  of  the  Ceram 
Sea  event  about  33  seconds  after  its  first  arrival. 
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Figure  IV- 2 displays  the  same  mixed  signal  data  in  the  1.5- 
2.  5 Hz  passband. 

B.  ADAPTIVE  BEAMFORMING  OUTPUTS 

The  composite  data  which  simulated  the  recorded  short-period 
data  at  KSRS  for  mixed  signal  situations  were  first  processed  by  the  adaptive 
beamformer  in  various  passbands.  To  display  the  wavetrains  for  the  signal 
from  the  Ceram  Sea,  we  processed  the  composite  data  by  beamforming  the 
array  using  the  exact  location  of  the  stronger  signal  source.  Time-varying 
adaptation  parameters  used  in  this  report  were  a=  0.  0045,  and  /3=  0.  65 
which  were  found  to  have  optimum  SNR  for  the  initial  P-wave.  Figure  IV- 3 
shows  the  beamforming  outputs  (0.  5-3.  5 Hz)  of  the  composite  data  by  steer- 
ing the  array  at  the  Ceram  Sea  event.  Because  the  signal  amplitude  from 
Ceram  was  much  larger  than  that  from  Tadzhik,  beamforming  of  the  stronger 
signal  would  eliminate  the  weaker  signal.  Therefore,  the  Tadzhik  signal  dis- 
appeared from  the  array  beams  in  this  case. 

To  demonstrate  the  weak  signal  extraction  by  suppressing  the 
stronger  interfering  signal  in  the  mixed- signal  situation.  Figure  IV- 4 shows 
the  beamforming  outputs  in  the  0.  5-1.  5 Hz  passband  of  the  composite  data 
where  the  mixed- signal  situation  is  simulated.  The  array  was  steered  at  the 
Tadzhik  location  (289°  and  14.  3 km/sec).  The  first  pulse  of  wavetrains  on 
the  ABF  output  is  the  initial  arrival  of  the  signal  from  the  Ceram  Sea.  The 
signal  from  Tadzhik  was  placed  33  seconds  after  the  initial  pulse  of  the  Ceram 
wavetrains  as  shown  on  the  ABF  beams.  The  figure  does  show  some  kind  of 
capability  that  the  adaptive  filter  removes  the  interfering  coda  and  success- 
fully separates  the  mixed- signals.  It  is  noted  that  a higher  fi  design  of  the 
ABF  would  totally  eliminate  the  Ceram  wavetrains,  but  it  also  degraded  the 
Tadzhik  signal  in  this  case. 
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STEERING  ARRAY  AT  TADZHIK 
(1.5-2.  5 Hz) 
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BEAMFORMING  OUTPUTS  OF  THE  MIXED- SIGNAL  DATA  BY 
STEERING  ARRAY  AT  TADZHIK,  USSR 
(0.5-3.  5 Hz) 
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where  fc  is  the  carrier  frequency  and  Af  is  one  half  of  a bandwidth.  Vari- 
ous bandwidths  of  data  were  tested  to  evaluate  the  effectiveness  of  single  fre- 
quency formulation. 

Figure  IV-7  illustrates  the  beamforming  outputs  of  the  mixed- 
signal  data  where  fc=  2.0  Hz,  and  Af  = 0.  3 Hz.  The  first  trace  is  the 
equally  weighted  array  beam  (beamsteer  output)  and  the  second  trace  is  the 
null  steered  array  beam  by  the  real -amplitude  weighting  algorithm.  The 
filter  weights  used  for  null  steering  was  listed  in  Figure  IH-3  where  the 
response  pattern  was  also  shown.  The  success  in  suppressing  the  'interfering' 
signed,  wavetrains  from  the  Ceram  Sea  to  noise  level  was  clearly  demonstrated 
in  the  figure. 

The  same  filter  weights  in  Figure  IH-3  were  applied  again  to 
the  same  mixed- signal  composte  data  with  different  bandwidths  in  bandpass 
filtering.  Figure  IV-8  shows  the  equally  weighted  array  beam  and  the  real- 
amplitude  weighted  array  beam  in  the  1.  0 Hz  bandwidth  (i.  e. , f = 2.  0 Hz, 
and  Af  = 0.  5 Hz).  In  this  case,  null  steering  was  still  very  effective  in  re- 
moving the  stronger  signal  wavetrains  from  the  Ceram  Sea  and  in  enhancing 
the  weaker  signal  from  Tadzhik.  Figure  IV -9  shows  the  relatively  wider 

band  data  with  f = 2.  0 Hz,  and  Af  = 1.  0 Hz.  Bandwidth  in  this  figure  was 
c 

doubled  from  that  of  Figure  IV-8,  and  effectiveness  and  validity  of  single 
frequency  formulation  was  degraded  somewhat  as  the  bandwidths  of  the  data 
increased.  However,  the  signal  from  Tadzhik  in  the  null  steered  array  beam 
in  Figure  IV-9  is  still  detectable,  while  in  the  equally  weighted  beam  it  is  not. 

The  complex  filter  weights  in  Table  IH-3,  obtained  by  complex- 

amplitude  shading  algorithm,  were  also  applied  to  the  same  sets  of  array  data 

for  various  equivalent  cases,  namely  f = 2.  0 Hz,  and  Af  = 0.  3,  0.  5,  and  1. 0 

c 

Hz.  Figure  IV-10  shows  the  complex- weighted  array  beams  where  the  first 
trace  is  the  beamsteer  output  (equally  weighted  with  zero  phase),  the  second 
trace  is  the  null- steered  beam  with  filter  weights  in  set  (a)  of  Table  IH-3; 
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and  the  third  trace  is  also  the  null-steered  beam  with  filter  weights  in  set  (b) 
of  the  same  table.  Both  filter  weight  constraint  conditions  had  the  same  re- 
sponse pattern  (Figure  III-5)  and  consistently  resulted  in  the  same  beam  out- 
puts (from  a visual  point  of  view).  Figure  IV -10  had  channel  data  input  the 
same  as  in  Figure  IV-7,  namely,  f = 2.0  Hz,  Af  = 0.  3 Hz.  However,  the 
complex  weighting  algorithm  in  this  case  failed  to  effectively  remove  the 
interfering  signal  from  the  Ceram  Sea  while  it  succeeded  in  the  real-amplitude 
algorithm  shown  in  Figure  IV-7.  However,  it  reduced  the  interfering  signal 
wavetrains  somewhat. 

The  same  complex  filter  weights  in  Table  III- 3 were  applied 
again  to  the  mixed-signal  composite  array  data  with  slightly  increased  band- 
width as  it  was  done  in  the  real -amplitude  algorithm.  Figure  IV-11  shows  the 
array  beams  in  which  bandpass  parameters  were:  f = 2.  0 Hz,  and  Af  = 0.  5 
Hz.  The  rejection  of  interfering  energy  deteriorated.  The  experiments  were 
carried  out  to  the  bandwidth  of  2.  0 Hz  signal  for  this  case  and  Figure  IV- 12 
displays  its  outputs  where  f = 2.  0 Hz  and  Af  =1.0  Hz.  The  experimental 
results  totally  failed  to  reduce  the  null- steered  signal  in  this  case. 

D.  SUMMARY 

Using  the  same  composite  set  of  the  stronger  interfering  sig- 
nal and  the  weaker  desired  signal  data,  the  experimental  work  of  mixed  sig- 
nal separation  was  carried  out  by  applying  the  ABF  processor  and  optimum 
array  shading  techniques  to  the  mixed- signal  data. 

The  adaptive  beamforming  of  time- domain  data  can  be  used  as 
a front-end  detector  or  estimator.  In  a detection  system,  only  azimuth  and 
velocity  are  needed  to  steer  the  array  which  can  be  fixed  to  cover  selected 
areas  of  interest.  Test  results  from  the  mixed  signal  case  indicated  that  the 
ABF  was  successful  in  eliminating  the  interfering  energy  and  in  extracting 
the  desired  weak  signal  to  the  level  that  a detection  can  be  claimed. 
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Experiments  in  various  passbands,  narrow  and  wide,  suggested  that  the  ABF 
processor  was  a powerful  method  for  a front-end  detection  system,  with  con- 
siderable capability  for  resolving  mixed  events. 

An  optimum  array  shading  technique  was  formulated  in  the  fre- 
quency domain  centered  about  a single  frequency.  The  method  can  be  used  as 
a post- detection  processor  where  the  velocity  and  azimuth  of  interfering  events 
are  known.  It  also  can  be  implemented  in  a detection  system  as  a special  pro- 
cessor to  reject  the  highly  coherent  noise.  Experimental  work  was  limited  to 
the  single-frequency  design  that  was  applied  to  a narrowband  signal.  The 
real-amplitude  algorithm  was  very  successful  in  removing  the  stronger  inter- 
fering signal  and  in  enhancing  the  desired  weaker  signal.  The  results  seem 
to  be  encouraging. 

However,  the  optimum  complex  weighting  technique  did  not 
yield  the  anticipated  results  for  which  it  was  designed.  Computation  in  this 
algorithm  is  laborious,  and  the  results  were  not  encouraging  at  all. 

In  short,  we  felt  that  the  ABF  processor  was  a powerful  beam- 
forming tool  with  the  capability  of  resolving  mixed  events.  Also,  the  real- 
amplitude  shading  algorithm  yielded  promising  results  in  our  experimental 
test.  The  experiments  with  the  ABF  and  shading  technique  for  separating 
mixed  events  are  considered  complimentary  in  that  the  shading  technique 
demonstrates  a physical  basis  for  the  performance  of  the  ABF  in  resolving 
mixed  events.  The  physical  basis  is  that  of  rapidly  time  varying  filter  weights 
applied  to  the  array  channel,  which  void  out  relatively  incoherent  energy  and 
pass  coherent  energy  in  accordance  with  Levin's  constraint. 
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SECTION  V 

CONCLUSION  AND  SUGGESTION 


A.  CONCLUSION 

The  physical  problem  of  mixed- signal  situations  we  have  tried 
to  solve  is  to  separate,  to  detect,  and  to  extract  the  weaker  signal  hidden  in 
the  stronger  wavetrains  of  an  earthquake.  Being  able  to  detect  weaker  sig- 
nals increases  the  detection  capability  of  an  array  and  being  able  to  extract 
lower  amplitude  signal  waveforms  from  much  higher  amplitude  interfering 
coherent  wavetrains  also  enhances  the  discrimination  capacity  of  the  array. 

We  have  developed  two  array  beamforming  methods  for  the 
mixed- signal  separation  problem.  One  is  the  time-domain  adaptive  beam- 
forming filter  and  the  other  the  frequency- domain  array  shading  (weighting) 
technique.  Both  methods  apply  the  variational  principle  of  mathematical 
physics.  However,  the  ABF  is  in  the  time-domain  and  data-iterative,  re- 
quires neither  a priori  noise  statistics  nor  signal  information  and  can  be 
implemented  in  real-time  processing.  The  array  shading  technique  is  in 
the  frequency-domain  design,  needs  information  of  interfering  azimuth,  and 
must  be  used  in  a post-detection  processor. 

Using  the  recorded  data  from  the  KSRS  short-period  array,  a 
composite  data  tape  was  formed  to  simulate  the  recorded  mixed  signal  data 
by  adding  the  data  of  high  amplitude  signal  wavetrains  from  a seismic  event 
in  the  Ceram  Sea  to  those  of  a low  amplitude  signal  (noise  level)  from  a 
seismic  event  in  Tadzhik.  The  ABF  processor  outputs  in  wide  band  (0.  5-3.  5 
Hz)  and  in  narrow  bands  (0.  5-1.  5 Hz,  and  1.  5-2.  5 Hz)  showed  that  the  inter- 
fering wavetrains  were  successfully  eliminated  and  the  desired  weak  signal 

V-l 


I 


I 

was  clearly  enhanced  to  a detectable  level,  while  the  desired  signal  in  beam- 
steer  outputs  was  at  the  level  of  interfering  wavetrains  and  was  non-detect- 
able.  The  experimental  test  of  the  real -amplitude  shading  algorithm  indicated 
that  the  results  were  promising  and  deserves  further  study.  However,  the 
complex-shading  algorithm  did  not  yield  the  satisfactory  results  predicted 
from  beam  pattern  computation.  Real  amplitude  weighting  appears  to  be  a 
sufficient  tool  for  studying  mixed- signal  problems. 

In  conclusion,  we  felt  that  the  present  ABF  is  a powerful  real- 
time array  processor  and  the  real -amplitude  shading  technique  could  be  use- 
ful for  post-processing  the  coda  of  large  signals. 


B.  SUGGESTION  FOR  FURTHER  STUDY 

The  ABF  could  be  improved  to  speed-up  the  computation. 

Technical  aspects  for  practical  implementation  as  a real-time  detector  are 
currently  being  studied.  A detailed  presentation  of  detection  performance  will 
be  in  a separate  report  (Shen,  1978)  where  it  was  tested  on  a large  ensemble 
of  events. 

The  preliminary  experimental  results  suggest  that  real- 
amplitude  shading  techniques  are  potentially  useful  for  removing  the  coherent 
interfering  energy  and  for  extracting  a hidden  signal  in  a coherent  interfer-  ^ 

ing  environment.  A useful  approach  would  be  to  bury  a signal  from  an  under- 
ground nuclear  explosion  in  various  types  of  interfering  environments.  A 
systematic  study  could  be  made  and  a set  of  filter  weights  could  be  designed 
from  the  study  so  that  they  could  be  used  for  post -detection  searches  for  hid- 
den explosions  or  other  research  purposes. 

The  real- amplitude  shading  algorithm  was  formulated  for  a 
single-frequency  signal  that  was  applied  to  a relatively  narrowband  signal. 


A wide  band  signal  can  be  synthesized  by  summing  the  data  of  various 


narrowband  signals  in  the  frequency  domain.  To  evaluate  its  practical  use- 
fulness, research  of  this  sort  should  be  conducted. 

Extension  of  the  present  work  to  long-period  seismic  array 
data  for  surface  wave  extraction  would  be  useful  for  both  detection  and  dis- 
crimination research. 
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